SUBROUTINE rouwenhorst(rho,sigma,p,y,f,n)
! Rowenhurst method to approximate univariate AR(1) process by Markov chain
!      y(t) = rho y(t-1)+ sigma sqrt(1-rho^2) e(t),   e(t)~N(0,1)
! INPUTS: rho - serial correlation coefficient,
!         sigma - coefficient of variation for y(t)
! OUTPUT: P is an n-by-n matrix of Markov transition probabilities
!         y is an n-by-1 vector of symmetric and evenly-spaced Markov state space
!         f is an n-by-1 vector of stationary distribution (binomial)
	IMPLICIT NONE

INTERFACE
SUBROUTINE grid(x,xmin,xmax,s,n)
    INTEGER, INTENT(IN) :: n,s
	DOUBLE PRECISION, DIMENSION(:), INTENT(OUT) :: x(n)
	DOUBLE PRECISION, INTENT(IN) :: xmin,xmax
	DOUBLE PRECISION :: c ! growth rate of grid subintervals for logarithmic spacing
	INTEGER :: i
END SUBROUTINE grid

SUBROUTINE binom(f,n)
	INTEGER, INTENT(IN) :: n
	DOUBLE PRECISION, DIMENSION(:), INTENT(OUT) :: f(n)
END SUBROUTINE binom
END INTERFACE


    INTEGER, INTENT(IN)    :: n
	DOUBLE PRECISION, INTENT(IN) :: rho, sigma
	DOUBLE PRECISION, DIMENSION(:,:), INTENT(OUT) :: p(n,n)
	DOUBLE PRECISION, DIMENSION(:), INTENT(OUT) :: y(n),f(n)
	DOUBLE PRECISION :: ybar, q
	IF (size(p,dim=1)/=n .or. size(p,dim=2)/=n) THEN
		PRINT '(a,i3,a,i3)', 'rouwenhorst: p must be a square matrix of size ',n,' x ',n
		STOP 'program terminated by rouwenhorst'
	END IF
	IF (size(f)/=n) THEN
		PRINT '(a,i3)', 'rouwenhorst: y and s must be vectors of the same size ',n
		STOP 'program terminated by rouwenhorst'
	END IF
	ybar=sigma*sqrt(real(n-1))
	q=(1+rho)/2
	CALL rhmat(p)
	CALL grid(y,-ybar,ybar,1,n)
	CALL binom(f,n)


CONTAINS
RECURSIVE SUBROUTINE rhmat(p)
	IMPLICIT NONE
	DOUBLE PRECISION, DIMENSION(:,:), INTENT(OUT) :: p
	DOUBLE PRECISION, DIMENSION(size(p,dim=1)-1,size(p,dim=2)-1) :: p1
	INTEGER :: h
	h=size(p,dim=1)
	IF (size(p,dim=2)/=h) STOP 'P must be a square matrix'
	IF (h<2) STOP 'P must be at least 2-by-2 matrix'
	IF (h==2) THEN
		p=reshape((/q,1-q,1-q,q/),(/2,2/))
	ELSE
		CALL rhmat(p1)
		p=0
		p(1:h-1,1:h-1)=q*p1
		p(1:h-1,2:h)=(1-q)*p1+p(1:h-1,2:h)
		p(2:h,1:h-1)=(1-q)*p1+p(2:h,1:h-1)
		p(2:h,2:h)=q*p1+p(2:h,2:h)
		p(2:h-1,:)=p(2:h-1,:)/2
	END IF
END SUBROUTINE rhmat
END SUBROUTINE rouwenhorst
